
#ldpc implementation
read "ldpc.mpl":

# compute the number of error-bits.
countErrorBits := proc(trans, corrected)
global N;
local i, result;
    result := 0;
    for i from 1 to N do
        if (trans[i] <> corrected[i]) then
            result := result + 1;
        end if;
    od;
    result
end proc:

# compute the bit error probability for 4-level cells
compErrorProb4 := proc(cells)
local cell, i, sum;
global N;
    sum := 0;
    for cell in cells do
        if (cell:-Bits[1] <> cell:-Bits0[1]) then
          sum := sum + 1;
	end if;
        if (cell:-Bits[2] <> cell:-Bits0[2]) then
          sum := sum + 1;
        end if;
    od;
    evalf(sum / N)
end proc;

# compute the cell level, e.g.,
# 00 -> 1, 01 -> 2, 10 -> 3, 11 -> 4
# FIXME: we currently do not use Gray codes
compCellLevel := proc(bits, bpc) option inline;
    1 + add(bits[i] * (2**(bpc - i)), i = 1..bpc)
end proc;

# compute the level of each cell from their initial bits stored
compCellLevels := proc(splitMsg, bpc) option inline;
    map(compCellLevel, splitMsg, bpc)
end proc;

# generates means of the log-resistances of all the cell levels.
genResistanceMeans := proc(start, gap, l) option inline;
    [seq(start+(i-1)*gap, i=1..l)]
end proc;

# compute a list of random variable following Gaussian distribution.
# each variable represents the distribution of the resistance of one
# cell level. The random variables in the returned list will be used
# for generating the random resistance of a cell.
genResistanceRandomVariableList := proc(l)
local means, dev, i;
    means := genResistanceMeans(3.0, 1.0, l);
    dev := 0.14;
    [seq(Statistics:-RandomVariable(Normal(means[i], dev)), i=1..l)]
end proc;

# the deviation of drift exponent = mean * SDMR
genDriftExponentRandomVariableList := proc(l)
local means, i, SDMR;
    SDMR := 0.4;
    # FIXME: right now we only support 4-level
    #means := [0.001, 0.02, 0.06, 0.10];
    #means := [0.005, 0.05, 0.08, 0.10];
    means := [0.05, 0.08, 0.09, 0.10];
    #means := [0.05, 0.08, 0.10, 0.12];
    [seq(Statistics:-RandomVariable(Normal(means[i], means[i]*SDMR)), i=1..l)]
end proc;

# counts the number of cells at each level from 1 to l
countCellsOfEachLevel := proc(cellLevels, l)
local results, i, cl;
    results := [seq(0, i = 1..l)];
    for cl in cellLevels do
        results := subsop(cl=(results[cl] + 1), results);
    od;
    results
end proc;

# Global random variables for generating random samples for resistances and
# drift exponents.
randvarsR := genResistanceRandomVariableList(4);
randvarsV := genDriftExponentRandomVariableList(4);

# computes the initial log-resistance and exponent of each cell
getInitialResistancesAndExponents := proc(cellLevels, l)
local cpl, cl; # cells per level
local resistanceSamples, exponentSamples; #[ [], [], ..., []];
local indexList, resistances, resistance, exponents, exponent, i;
global randvarsR, randvarsV;

    cpl := countCellsOfEachLevel(cellLevels, l);
    resistanceSamples := [seq(convert(Statistics:-Sample(randvarsR[i], cpl[i]), list), i=1..l)];
    exponentSamples := [seq(convert(Statistics:-Sample(randvarsV[i], cpl[i]), list), i=1..l)];
    # distribute resistance to each cell
    indexList := [seq(1, i=1..l)];
    resistances := [];
    exponents := [];
    for cl in cellLevels do
        i := indexList[cl];
        resistance := (resistanceSamples[cl])[i];
        resistances := [op(resistances), resistance];
        exponent := (exponentSamples[cl])[i];
        exponents := [op(exponents), exponent];
        indexList[cl] := i + 1;
    od;
    [resistances, exponents]
end proc;

# generates a new sample of initial resistance.
getNewResistance := proc(l) option inline;
global randvarsR;
    (convert(Statistics:-Sample(randvarsR[l], 1), list))[1]
end proc;

# generates a new sample of drift exponent.
getNewDriftExponent := proc(l) option inline;
global randvarsV;
    (convert(Statistics:-Sample(randvarsV[l], 1), list))[1]
end proc;

# given a list of bits to store, we generate a list of records, where
# each record represents a cell, and each record field represents a
# property associated with the cell. A cell currently has the
# following properties:
# 0. initial log-resistant (log base is 10)
# 1. current log-resistant
# 2. current drifting exponent
# 3. time elapsed since the cell was written
# 4. the bits which are represented by the current resistant
# 5. the levels which correspons to the bits
initPCMCells := proc(msg, l)
local size, bpc; #bits per cell
local cellnum, splitMsg;
local level, cellLevels, re;
local i, cell, rs, es;
    size := nops(msg);
    bpc := sqrt(l);
    if (whattype(bpc) <> integer) then
        error "User Error: l is required to be power of 2";
    end if;
    if (irem(size, bpc) <> 0) then
        error "User Error: the size of message bits can not be divided exactly by bpc";
    end if;
    cellnum := iquo(size, bpc);
    splitMsg := [ListTools:-LengthSplit(msg, bpc)];
    cellLevels := compCellLevels(splitMsg, bpc);
    re := getInitialResistancesAndExponents(cellLevels, l);
    rs := re[1];
    es := re[2];
    # for i from 1 to cellnum do
    #     cell := Record('R0' = rs[i], 'R' = rs[i], 'V' = es[i], 'T' = 0, 'Bits' = splitMsg[i], 'L' = cellLevels[i],  'Bits0' = copy(splitMsg[i]), 'L0' = cellLevels[i]);
    #     cells := [op(cells), cell];
    # od;
    [seq(Record('R0' = rs[i], 'R' = rs[i], 'V' = es[i], 'T' = 0, 'Bits' = splitMsg[i], 'L' = cellLevels[i],  'Bits0' = copy(splitMsg[i]), 'L0' = cellLevels[i]), i = 1..cellnum)]
end proc;

# plot the histograms of the log-resistances and drift exponents of the
# cells with level l
plotOneLevel := proc(cells, l)
local rs, sdev, cell; # for resistances
local vs, sdevV; #for drift exponents
    rs := [];
    vs := [];
    for cell in cells do
        if (cell:-L = l) then
            rs := [op(rs), cell:-R];
            vs := [op(vs), cell:-V];
        end if;
    od;
    sdev := Statistics:-StandardDeviation(rs);
    sdevV := Statistics:-StandardDeviation(vs);
    <Statistics:-Histogram(rs, binwidth=sdev, frequencyscale = absolute), Statistics:-Histogram(vs, binwidth=sdevV, frequencyscale = absolute, color = red)>;
#    <Statistics:-Histogram(rs, frequencyscale = absolute), Statistics:-Histogram(vs, frequencyscale = absolute, color = red)>;
end proc:

# plot the histograms of the log-resistances and drift exponents of the
# cells with levels from 1 to l
plotAll := proc(cells, l)
local ps, i;
    ps := <plotOneLevel(cells, 1)>;
    for i from 2 to l do
        ps := < ps | plotOneLevel(cells, i) >;
    od;
    plots:-display(ps)
end proc:

# compute the threshold between level i and i+1.
compThreshold := proc(i, means, devs) option inline;
    (means[i+1]*devs[i] + means[i] * devs[i+1])/(devs[i] + devs[i+1])
end proc;

# compute the initial thresholds for each two adjacent levels
compThresholds := proc(l)
local means, devs;
    means := [3.0, 4.0, 5.0, 6.0]; # genResistanceMeans(3.0, 1.0, l);
    devs := [seq(0.14, i=1..l)];
    [seq(compThreshold(i, means, devs), i=1..(l-1))]
end proc;

# thresholds for 4-level PCM cells
cellLevelThresholds := compThresholds(4):

# the log-based version of the cell drift model.
# returns the drifted resistance only
driftOneCell := proc(cell) option inline;
    evalf(cell:-R0 + (cell:-V) * log[10](cell:-T))
end proc:

# compute the new level given a new resistance and a pre-determined
# threshold set.
compNewLevel := proc(r, thresholds)
local level, threshold;
    level := 1;
    for threshold in thresholds do
        if r < threshold then
            return level;
        end if;
        level := level + 1;
    od;
    level
end proc:

# obtain the corresponding bit list from the cell level.
level2Bits4 := proc(level)
    if level = 1 then
        return [0, 0];
    elif level = 2 then
        return [0, 1];
    elif level = 3 then
        return [1, 0];
    else
        return [1, 1];
    end if;
end proc;

# drift all the cells
# via 1. updating resistance R, 2. update the elapsed time of each
# cell, 3. update the corresponding bits and the levels
# FIXME: we assume that we are using 4-level cell
driftAllCells := proc(cells, t)
local cell, r, newlevel;
global cellLevelThresholds;
    # assume we are using 4-level cells
    for cell in cells do
        cell:-T := cell:-T + t;
        cell:-R := driftOneCell(cell);
        newlevel := compNewLevel(cell:-R, cellLevelThresholds);
        if (cell:-L <> newlevel) then
            cell:-L := newlevel;
            cell:-Bits := level2Bits4(newlevel);
        end if;
    od;
end proc:

# drift all the cells
# via 1. updating resistance R, 2. update the elapsed time of each
# cell, 3. update the corresponding bits and the levels
# FIXME: we assume that we are using 4-level cell
justComputeNewLevels := proc(cells)
local cell, r, newlevel;
global cellLevelThresholds;
    # assume we are using 4-level cells
    for cell in cells do
        newlevel := compNewLevel(cell:-R, cellLevelThresholds);
        if (cell:-L <> newlevel) then
            cell:-L := newlevel;
            cell:-Bits := level2Bits4(newlevel);
        end if;
    od;
end proc:




# put bit list in each cell back into one list
collectBitsIntoOne := proc(cells) option inline;
    [seq(op(c:-Bits), c = cells)]
end proc;

# compute the level whose corresponding bits' i-th bit has value v
# FIXME: only 4 level is supported
getDesiredLevels := proc(i, v)
local levels, j, result;
    levels := [[0, 0], [0, 1], [1, 0], [1, 1]];
    result := [];
    for j from 1 to 4 do
        if (levels[j])[i] = v then
            result := [op(result), j];
        end if;
    od;
    result
end proc;



# cache for storing means and devs for R at level i time t.
rdCache := table();
rdkeys := {};

# computes mu_i(t) and (delta_i^2)(t)
# i: cell level
# t: time
# l: number of cell levels
# the function returns a pair (mu, dev^2)
compNewRD := proc(i, t, l)
local meansR, meansV, devR, devV, SDMR, lgt, result;
global rdCache, rdkeys;
    if member([i, t], rdkeys) then
        return  rdCache[(i, t)];
    else
        SDMR := 0.4;
        meansR := [3.0, 4.0, 5.0, 6.0]; #genResistanceMeans(3.0, 1.0, l);
        # FIXME: right now we only support 4-level
        #meansV := [0.001, 0.02, 0.06, 0.10];
        #meansV := [0.005, 0.05, 0.08, 0.10];
        #meansV := [0.05, 0.06, 0.09, 0.10];
        #meansV := [0.05, 0.08, 0.10, 0.12];
        meansV := [0.05, 0.08, 0.09, 0.10];
        devR := 0.14;
        devV := meansV[i] * SDMR;
        lgt := evalf(log[10](t));
        result := [meansR[i] + meansV[i] * lgt, devR^2 + (devV^2) * (lgt^2)];
        rdCache[(i, t)] := result;
        rdkeys := rdkeys union {[i, t]};
        return result;
    end if;
end proc;

# buggy, don't use for now. Only for debugging.
# compTi := proc(i, t)
#     rdi := compNewRD(i, t, 4);
#     rdi1 := compNewRD(i+1, t, 4);
#     evalf((rdi1[1] * rdi[2] + rdi[1] * rdi1[2])/(rdi[2] + rdi1[2]))
# end proc;


# manually evaluate PDF
evalPDF := proc(md, r) option inline;
    exp(((r - md[1])^2)/(-2*md[2]))/sqrt(2*Pi*md[2])
end proc:

# computePThreadSafe := proc(r, t, l, i, k)
# local levels;
#     levels := getDesiredLevels(i, k);
#     # use unrolled version.
#     (evalPDF(compNewRD(levels[1], t, l), r) + evalPDF(compNewRD(levels[2], t, l), r))/l
#    #add(evalPDF(compNewRD(level, t, l), r), level = getDesiredLevels(i, k))/l
# end proc:
computePThreadSafe := proc(r, t, l, levels) option inline;
    # use unrolled version.
    (evalPDF(compNewRD(levels[1], t, l), r) + evalPDF(compNewRD(levels[2], t, l), r))/l
   #add(evalPDF(compNewRD(level, t, l), r), level = getDesiredLevels(i, k))/l
end proc:

# a faster implementation of computeP0P1slow
computeP0P1unroll := proc(r, t, l)
local cp01, cp11, cp02, cp12;
    cp01 := computePThreadSafe(r, t, l, [1, 2]);
    cp11 := computePThreadSafe(r, t, l, [3, 4]);
    cp02 := computePThreadSafe(r, t, l, [1, 3]);
    cp12 := computePThreadSafe(r, t, l, [2, 4]);
    # cp01 := computePThreadSafe(r, t, l, 1, 0);
    # cp11 := computePThreadSafe(r, t, l, 1, 1);
    # cp02 := computePThreadSafe(r, t, l, 2, 0);
    # cp12 := computePThreadSafe(r, t, l, 2, 1);
    evalf(cp11/(cp01 + cp11)), evalf(cp12/(cp02 + cp12))
end proc:

# compute probabilities of variable nodes for the bits stored in one
# cell.
compInitialProb := proc(cell, l)
local t;
    t := cell:-T;
    if (t = 0) then
        return op(cell:-Bits);
    else
        return computeP0P1unroll(cell:-R, t, l)
    end if;
    #computeP0P1(cell:-R, t, l, 1), computeP0P1(cell:-R, t, l, 2)
end proc;

# a slightly faster function for computing initial probabilities of
# variable nodes. We use the experimental Maple task model to
# parallelizing sequence construction
compInitialProbs := proc(cells, l) option inline;
    #Array([seq(compInitialProb(i, l), i = cells)])
    #[seq(compInitialProb(i, l), i = cells)]
    map(compInitialProb, cells, 4)
end proc;

# run ECC to find out errors.
errorCorrect := proc(cells, l, H, itermax, tGraph)
local msg, vnodeProbs; # initial probabilities for variable nodes.
    msg := collectBitsIntoOne(cells);
    vnodeProbs := compInitialProbs(cells, l);
    myprint("BER before applying EC: ", compErrorProb4(cells));
    decodeSPA(msg, H, itermax, vnodeProbs, tGraph);
end proc:

# periodically scrub the memory cells
# msg : codewords to be stored
# l : cell level
# cps: number of cells in each segment
# gap: time gap between scrubbing of two different segments
# maxtime : time to stop measuring.
# itermax : iteration limit for LDPC decoding
# H: parity check matrix for LDPC code
# opt: switch to turn on optimizaiton to utilize error log before
# applying ECC.
periodicalScrub := proc(msg, l, cps, gap, maxtime, itermax, H, opt)
local tGraph; # for storing the Tanner graph created from H.
local t; # time
local s; # number of segments
local numcell;
local i; # index of the segment that is being scrubbed
local j, k, cindex, bindex, cellBits, re;
local errorLog; # a set of pair [index, the faulty bit value]
local cells, bpc;
local currentBits, correctedBits, cellsToRefresh; # for storing the corrected bits given by ECC, and the bits before correction.
# the following locla variables are performance counters.
local bers; # a list of pair stores the bit error rate of the block before and after running ECC.
local correctblkno, undetblkno, detblkno; # number of correct, undetected, detected blocks.
local totalErrorBlocks, totalIters, successflag;
local berbefore, berafter;
local decodedResult;
local nosegscrubbed, stat;
local nofaultybits;
local bitsBforeECC;
local baddataname;
global savedir;

    tGraph := TannerGraph(H);

    # initializing performance counters.
    bers := [];
    correctblkno := 0;
    undetblkno := 0;
    detblkno := 0;
    totalErrorBlocks := 0;
    totalIters := 0;
    nosegscrubbed := 0;
    nofaultybits := 0;



    cells := initPCMCells(msg, l);
    bpc := nops(cells[1]:-Bits);
    errorLog := {};
    numcell := nops(cells);
    t := gap;  # after writing the whole codeword for gap time, we start scrubbing
    driftAllCells(cells, gap);
    i := 1;
    s := iquo(numcell, cps);


    if opt then
        print("Use error log to correct errors before ECC");
    else
        print("Don't use error log before ECC");
    end if;

    myprint("Total number of segments: ", s);
    myprint("Cells per segment: ", cps);

    while t <= maxtime do
        myprint("Segment being scrubbed: ", i);
        myprint("Current time: ", t);

        if opt then
            # step 1.1 correct faulty bits found in the previous round
            cellsToRefresh := {}; # absolute indices of cells to be refreshed before performing LDPC decoding
            bitsBforeECC := 0;
            if nops(errorLog) > 0 then
                for j from 1 to cps do
                    cindex := (i - 1) * cps + j;
                    cellBits := cells[cindex]:-Bits;
                    for k from 1 to bpc do
                        bindex := (cindex - 1) * bpc + k;
                        # if the bit is still faulty i.e., didn't change back
                        # to correct bit, then start correcting.
                        if member([bindex, cellBits[k]], errorLog) then
                            bitsBforeECC :=  bitsBforeECC + 1;
                            # correct the faulty bit
                            (cells[cindex]:-Bits)[k] := (cellBits[k] + 1) mod 2;
                            # keep in mind that we need to refresh this cell
                            # given the newly stored bits.
                            cellsToRefresh := cellsToRefresh union {cindex};
                        end if;
                    od;
                od;
            end if;
            print("Number of faulty bits corrected before applying EC: ", bitsBforeECC);
            print("Total number of faulty bits in error log before applying EC: ", nops(errorLog));

            # step 1.2 refresh the cells which store the faulty bits
            for cindex in cellsToRefresh do
                # update the initial resistance, current resistance,
                # exponent from the newly computed level.  The reset
                # elapsed time to 0.
                cells[cindex]:-L := compCellLevel(cells[cindex]:-Bits, bpc);
                cells[cindex]:-T := 0;
                cells[cindex]:-R0 := getNewResistance(cells[cindex]:-L);
                cells[cindex]:-R := cells[cindex]:-R0;
                cells[cindex]:-V := getNewDriftExponent(cells[cindex]:-L);
            od;
        end if;


        # save the ber before ECC.
        berbefore := compErrorProb4(cells);

        # step 2. run ECC (LDPC here)
        decodedResult := errorCorrect(cells, l, H, itermax, tGraph);
        successflag := decodedResult[1];
        correctedBits := decodedResult[3];

        # if decoding doesn't terminate, take a snapshot of the cells
        # and the probabilities.
        if not successflag then
            baddataname:= cat(savedir, cat( cat(StringTools:-FormatTime("%m-%d-%H-%M-%S"), ".snap.mpl")));
            save cells, baddataname;
        end if;

        # if the decoding is successful we refresh
        # or if the decodig is not terminating, and the error log is
        # empty i.e., there's no more good hope, we use best guess.
        if (successflag) or ((not successflag) and nops(errorLog) = 0) then
            if opt then
                # cleanup the error log.
                errorLog := {};
                currentBits := collectBitsIntoOne(cells);
                for j from 1 to nops(msg) do
                    if correctedBits[j] <> currentBits[j] then
                        errorLog := errorLog union {[j, currentBits[j]]};
                        #errorLog := {op(errorLog), [j, currentBits[j]]};
                    end if;
                od;
            end if;

            # step 3. If LDPC terminates, update errorLog with the new
            # error information clean up all the cells in the current
            # segment by writing the bits returned by LDPC decoding
            # back to the segment
            for j from 1 to cps do
                cindex := (i - 1) * cps + j;
                cellBits := copy(cells[cindex]:-Bits);
                # correct all the bits in the cell cindex
                for k from 1 to bpc do
                    bindex := (cindex - 1) * bpc + k;
                    # correct the faulty bit
                    (cells[cindex]:-Bits)[k] := correctedBits[bindex];
                    # remove the entry from errorLog if this bit is found
                    # faulty by ECC
                    if opt then
                        if member([bindex, cellBits[k]], errorLog) then
                            errorLog := errorLog minus {[bindex, cellBits[k]]};
                        end if;
                    end if;
                od;
                # refresh this cell
                # update the resistance, exponent from the newly computed level.
                # The reset elapsed time to 0.
                cells[cindex]:-L := compCellLevel(cells[cindex]:-Bits, bpc);
                cells[cindex]:-T := 0;
                cells[cindex]:-R0 := getNewResistance(cells[cindex]:-L);
                cells[cindex]:-R := cells[cindex]:-R0;
                cells[cindex]:-V := getNewDriftExponent(cells[cindex]:-L);
            od;
        else
            print("Decoding result is not reliable, while error log still has correct bits, we don't do anything after this scrubbing");
        end if;

        berafter := compErrorProb4(cells);
        myprint("BER after applying ECC: ", berafter);

        bers := [op(bers), [berbefore, berafter]];
        totalIters := totalIters + decodedResult[2];
        nosegscrubbed := nosegscrubbed + 1;
        nofaultybits := countErrorBits(msg, correctedBits);
        if (nofaultybits = 0) then
            if successflag then
                correctblkno := correctblkno + 1;
            else # this branch is not quite possible to happen unless a bug.
                error "warning, a bug happens";
            end if;
        else #if errbits > 0 then
            totalErrorBlocks := totalErrorBlocks + 1;
            if successflag then
                # then this is an undetected error
                undetblkno := undetblkno + 1;
            else
                detblkno := detblkno + 1;
            end if;
        end if;

        # increase time and move to the next segment (by increasing
        # the index i)
        # if this is the last segment, go back to scrub the first one
        if i = s then
            i := 1;
        else
            i := i + 1;
        end if;
        t := t + gap;
        # drift all the cells according to time gap
        driftAllCells(cells, gap);
    od;
    print("Statistics of periodical scrubbing");
    stat := [evalf(totalIters/nosegscrubbed), nosegscrubbed, correctblkno, detblkno, undetblkno, nofaultybits, bers];
    print(stat);
    #collectBitsIntoOne(cells);
    stat
end proc:

# without any scrubbing.
noScrub := proc(msg, l, maxtime, itermax, H)
local cells, tGraph;
local recvmsg, vnodeProbs; # initial probabilities for variable nodes.
    cells := initPCMCells(msg, l);
    tGraph := TannerGraph(H);
    driftAllCells(cells, maxtime);
    recvmsg := collectBitsIntoOne(cells);
    vnodeProbs := compInitialProbs(cells, l);
    myprint("BER before applying ECC: ", compErrorProb4(cells));
    decodeSPA(recvmsg, H, itermax, vnodeProbs, tGraph);
end proc:

# the naive scrubbing. That is, the scrubbing is periodical, but it
# cleansup all the bits at the same time. If the scrubbing does not
# terminate then do not cleanup any bits.
basicScrub := proc(msg, l, gap, maxtime, itermax, H)
local tGraph; # for storing the Tanner graph created from H.
local t; # time
local numcell; # total number of cells
local i; # index of the segment that is being scrubbed
local j, k, cindex, bindex, cellBits, re;
local cells, bpc;
local currentBits, correctedBits, cellsToRefresh; # for storing the corrected bits given by ECC, and the bits before correction.
# the following locla variables are performance counters.
local bers; # a list of pair stores the bit error rate of the block before and after running ECC.
local correctblkno, undetblkno, detblkno; # number of correct, undetected, detected blocks.
local totalErrorBlocks, totalIters, successflag;
local berbefore, berafter;
local decodedResult;
local scrubbedtimes, stat;
local nofaultybits;
local bitsBforeECC;

    tGraph := TannerGraph(H);

    # initializing performance counters.
    bers := [];
    correctblkno := 0;
    undetblkno := 0;
    detblkno := 0;
    totalErrorBlocks := 0;
    totalIters := 0;
    scrubbedtimes := 0;
    nofaultybits := 0;

    cells := initPCMCells(msg, l);
    bpc := log[2](l);
    numcell := nops(cells);
    t := gap;  # after writing the whole codeword for gap time, we start scrubbing
    driftAllCells(cells, gap);

    while t <= maxtime do
        myprint("Current time: ", t);

        # save the ber before ECC.
        berbefore := compErrorProb4(cells);

        # step 2. run ECC (LDPC here)
        decodedResult := errorCorrect(cells, l, H, itermax, tGraph);
        successflag := decodedResult[1];
        correctedBits := decodedResult[3];
#        if successflag then
        # step 3. If LDPC terminates, clean up all the cells by
        # writing the bits returned by LDPC decoding back to the
        # segment
        # we clean up the segment no matter what.
        for cindex from 1 to numcell do
            # correct all the bits in the cell cindex
            for k from 1 to bpc do
                bindex := (cindex - 1) * bpc + k;
                # correct the faulty bit
                (cells[cindex]:-Bits)[k] := correctedBits[bindex];
            od;
            # refresh this cell
            # update the resistance, exponent from the newly computed level.
            # The reset elapsed time to 0.
            cells[cindex]:-L := compCellLevel(cells[cindex]:-Bits, bpc);
            cells[cindex]:-T := 0;
            cells[cindex]:-R0 := getNewResistance(cells[cindex]:-L);
            cells[cindex]:-R := cells[cindex]:-R0;
            cells[cindex]:-V := getNewDriftExponent(cells[cindex]:-L);
        od;
        # else
        #     print("Decoding result is not reliable, we don't do anything after this scrubbing");
        # end if;

        berafter := compErrorProb4(cells);
        myprint("BER after applying ECC: ", berafter);

        bers := [op(bers), [berbefore, berafter]];
        totalIters := totalIters + decodedResult[2];
        scrubbedtimes := scrubbedtimes + 1;
        nofaultybits := countErrorBits(msg, correctedBits);
        if (nofaultybits = 0) then
            if successflag then
                correctblkno := correctblkno + 1;
            else # this branch is not quite possible to happen unless a bug.
                error "warning, a bug happens";
            end if;
        else #if errbits > 0 then
            totalErrorBlocks := totalErrorBlocks + 1;
            if successflag then
                # then this is an undetected error
                undetblkno := undetblkno + 1;
            else
                detblkno := detblkno + 1;
            end if;
        end if;

        # increase time
        t := t + gap;
        # drift all the cells according to time gap
        driftAllCells(cells, gap);
    od;
    print("Statistics of periodical scrubbing");
    stat := [evalf(totalIters/scrubbedtimes), scrubbedtimes, correctblkno, detblkno, undetblkno, nofaultybits, bers];
    print(stat);
    #collectBitsIntoOne(cells);
    stat
end proc:



################### performance measurement functions #######################
# count the number of erroneous bits whose resistance level has become
# larger than the threshold of that level.
countErrors := proc(cells, l)
local thresholds, cell, level, counters;
    thresholds := compThresholds(l);
    counters := [seq(0, i=1..(l-1))];
    for cell in cells do
        level := cell:-L0;
        if level <> l then
            if cell:-R > thresholds[level] then
                counters[level] := counters[level] + 1;
            end if;
        end if;
    od;
    counters
end proc;



# find the probability which is not too large, i.e. low < p < high
findUsefulPoints := proc(probs, high, low)
local results, prob;
    results := [];
    for prob in probs do
        if prob < high and prob > low then
            results := [op(results), prob];
        end if;
    od;
    results
end proc;

# a test for initPCMCells
msg := [seq(rand() mod 2, i = 1..4376)]:
l := 4;
cells := initPCMCells(msg, l):
#plotAll(cells, l);


# test drifting in sec
# 10**7 = 115 days
#driftTime := (10**6);
#driftAllCells(cells, driftTime):
# plotAll(cells, l):
# countErrors(cells, l);
#compErrorProb4(cells);
# probs := compInitialProbs(cells, l):
# nops(findUsefulPoints(probs, 0.8, 0.2));

# compute raw error probability at time point n * gap
rawErrorProbs := proc(gap, n)
local probs, msg, l, i, cells;
    msg := [seq(rand() mod 2, i = 1..1998)]:
    l := 4;
    cells := initPCMCells(msg, l):
    probs:= [];
    for i from 1 to n do
        driftAllCells(cells, gap);
        probs := [op(probs), compErrorProb4(cells)];
    od;
    probs;
end proc:

# compute average raw error probability at time point n * gap
plotAverageRowEBRs := proc(k, gap, n)
local probs, points, i;
    probs := [seq(0, i=1..n)];
    for i from 1 to k do
        probs := probs + rawErrorProbs(gap, n);
    od;

    points := [];
    for i from 1 to n do
        points := [op(points), [i, evalf(probs[i]/k)]];
    od;
    plot(points, style=line)
end proc:


getCellsOfSegment := proc(cells, s, cps)
    [seq(cells[i], i = ((s-1) * cps + 1)..(s * cps))]
end proc:

# count the number of 0s and 1s stored in a list of PCM cells.
countZeroOne := proc(cells)
local counters, bits, b;
    counters := [0,0];
    bits := collectBitsIntoOne(cells);
    for b in bits do
        if b = 1 then
            counters[1] := counters[1] + 1;
        else
            counters[2] := counters[2] + 1;
        end if;
    od;
    counters
end proc:

# count the number of zeros and ones in a list of cells with
# s segments and cps cells per segment.
countAllZeroOne := proc(cells, s, cps)
    [seq(countZeroOne(getCellsOfSegment(cells, i, cps)), i = 1..s)]
end proc:
